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Abstract We present recent results on coarse-graining techniques for thermody- 
namic quantities (canonical averages) and dynamical quantities (averages of path 
functionals over solutions of overdamped Langevin equations). The question is how 
to obtain reduced models to compute such quantities, in the specific case when the 
functional to be averaged only depends on a few degrees of freedom. We mainly 
review, numerically illustrate and extend results from |[3][18), concerning the com- 
putation of the stress-strain relation for one-dimensional chains of atoms, and the 
construction of an effective dynamics for a scalar coarse-grained variable when the 
complete system evolves according to the overdamped Langevin equation. 



1 Motivation 

In molecular simulation, two types of quantities are typically of interest: averages 
with respect to the canonical ensemble (thermodynamic quantities, such as stress, 
root-mean-square distance, . . . ), and averages of functionals over paths (dynamic 
quantities, like viscosity, diffusion coefficients or rate constants). In both cases, the 
question of coarse-graining is relevant, in the sense that the considered function- 
als typically depend only on a few variables of the system (collective variables, or 
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reaction coordinates). Therefore, it is essential to understand how to obtain coarse- 
grained models on these variables. 



1.1 Coarse-graining of thermodynamic quantities 

Computing canonical averages is a standard task in molecular dynamics. For a 
molecular system whose atom positions are described by a vector q € M", these 
quantities read 

/ &(q)dll 

where <t> : R" — > R is the observable of interest and /i is the Boltzmann-Gibbs mea- 
sure, 

d\i =Z- 1 exp{-pV(q))dq, (1) 

where V is the potential energy of the system, j3 is proportional to the inverse of 

the system temperature, and Z= / exp(— j5V(q))dq is a normalizing constant. 

Typically, q represents the position of N particles in dimension d, hence q G R" 
with n = dN. 

As mentioned above, observables of interest are often functions of only part of 
the variable q. For example, q denotes the positions of all the atoms of a protein 
and of the solvent molecules around, and the quantity of interest is only a particu- 
lar angle between some atoms in the protein, because this angle characterizes the 
conformation of the protein (and thus the potential energy well in which the system 
is, is completely determined by the knowledge of this quantity of interest). Another 
example is the case when q = (q l ,. .. ,q") denotes the positions of all the atoms of 
a one-dimensional chain, and quantities of interest are only a function of the total 
length q" — q l of the chain. 

We thus introduce the so-called reaction coordinate 

% : R" -> R, 

which contains all the information we are interested in. Throughout this article, we 
assume that it is a smooth function such that |V|| is bounded from below by a 
positive constant, so that the configurational space can be foliated by isosurfaces 
associated to t, . A simple case that will be considered below is t, (q l , . . . , q") = q". 

To this function ^ is naturally associated an effective energy A, called the free 
energy, such that 

d{% *n) = exp(-j3A(z))<fe, 

where £, * ji denotes the image of the measure jU by % . In other words, for any test 
function 4> : R ->• R, 



/ <P($(q))Z- l exp(-l3V(q))dq= f <P(z) exp(- j3A(zj) dz. (2) 
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Expressions of A and its derivative are given below (see Section [T~4l i. 

The interpretation of (O is that, when Q is a random variable distributed accord- 
ing to the Boltzmann measure (Q}, then | (Q) is distributed according to the measure 
exp(— j3A (z)) dz. Hence, the free energy A is a relevant quantity for computing ther- 
modynamic quantities, namely canonical averages. 

In conclusion, the question of coarse-graining thermodynamic quantities amounts 
to computing the free energy, and there are several efficient methods to perform such 
calculations (see for example J6][l9l). In the sequel of this article, we address a par- 
ticular case, motivated by materials science, where the system under consideration 
is a one-dimensional chain of atoms, and %(q l , . . . ,q n ) = q" — q l is the length of 
the chain (see Fig. Q] below). We are interested in the free energy associated to this 
reaction coordinate, and its behaviour when the number n of particles goes to +°°. 
Standard algorithms to compute the free energy then become prohibitively expen- 
sive, as the dimension of the system becomes larger and larger. Alternative strategies 
are needed, and we investigate analytical methods, based on large deviations princi- 
ples, in Section[2] 



1.2 Coarse-graining of dynamical quantities 

The second topic of this contribution is related to the dynamics of the system, and 
how to coarse-grain it. In short, we will show how to design a dynamics that ap- 
proximates the path t\-^% (Qt), where £, is the above reaction coordinate. 

To make this question precise, we first have to choose the full dynamics, which 
will be the reference one. In the following, we consider the overdamped Langevin 
dynamics on state space W: 



where W t is a standard 77 -dimensional Brownian motion. Under suitable assumptions 
on V, this dynamics is ergodic with respect to the Boltzmann-Gibbs measure (Q]i 
(see J3 and references therein). Hence, for /i-almost all initial conditions Qq, 



almost surely. In practice, this convergence is often very slow, due to some metasta- 
bilities in the dynamics: Q, samples a given well of the potential energy for a long 
time, before hoping to some other well of V. 

An important dynamical quantity we will consider below is the average residence 
time, that is the mean time that the system spends in a given well, before hoping to 
another one, when it follows the dynamics (01. Typically, the wells are fully de- 
scribed through 4 {q is in a given well if and only if £, (q) is in a given interval), so 
that these times can be obtained from the knowledge of the time evolution of £, (Qt), 
which is expensive to compute since it means simulating the full system. 



dQ t = -VV{Q t )dt + ^/lf ri dW t , Q t=0 = Q 



(3) 




(4) 
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In Section|3]below, we will first present a one-dimensional dynamics of the form 
drj t = b{n t )dt + ^- 1 a%)dB„ (5) 

where B, is a standard one-dimensional Brownian motion and b and a are scalar 
functions, such that (f[ t ) 0<t<T i s a good approximation (in a sense to be made 
precise below) of (£, (2r))o<r<r- Hence, the dynamics © can be thought of as a 
coarse-grained, or effective, dynamics for the quantity of interest. A natural require- 
ment is that © preserves equilibrium quantities, i.e. it is ergodic with respect to 
exp(— j5A(z))dz, the equilibrium measure of £,(Q t ) when Q, satisfies (f3]), but we 
typically ask for more than that. For example, we would like to be able to recover 
residence times in the wells from (©, hence bypassing the expensive simulation of 

As a matter of fact, the coarse-grained dynamics 

dz t = -A'{z,)dt + dB„ (6) 

is a one-dimensional dynamics that is ergodic with respect to exp(— fiA(z))dz. It 
can thus be thought of as a natural candidate for a dynamics approximating ^ (Qt), 
all the more so as practitioners often look at the free energy profile (i.e. the func- 
tion z i-> A(z)) to get an idea of the dynamics of transition (typically the transi- 
tion time) between one region indexed by the reaction coordinate (say for example 
{q € K"; 4 (q) < zo}) and another one (for example {q G K"; h, (q) > zo}). If % (Qt) 
follows a dynamics which is close to @, then the Transition State Theory says that 
residence times are a function of the free energy barriers |fT7l[T6l , and then it makes 
sense to look at the free energy to compute some dynamical properties. It is thus 
often assumed that there is some dynamical information in the free energy A. 

In the sequel, we will compare the accuracy of both coarse-grained dynamics, an 
effective dynamics of type (© (namely dynamics d65l l below) and the dynamics (01 
driven by the free energy. Their relation has been investigated from an analytical 
viewpoint in (H Section 2.3] (see also 03] Sec. 10 and Eq. (89)] and |20l ). 



1.3 Outline of the article 

In this contribution, we mainly review, numerically illustrate and extend results from 
the two articles |3][T8). The aim is to present in a pedagogical and unified manner 
recent contributions on coarse-graining procedures concerning: (i) a static case in- 
spired by material sciences, namely the computation of stress-strain relation for 
one-dimensional chains of atoms, in the thermodynamic limit (Section [U and (ii) 
a dynamic case inspired by molecular dynamics computations, namely the deriva- 
tion of effective dynamics along the reaction coordinate, for overdamped Langevin 
equations (Section|3]l. Compared to the original articles j3][T8l, we propose some ex- 
tensions of the theoretical results (see e.g. Section lZ2l ). some simpler proofs in more 
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restricted settings (in Section 13 . 31 > and new numerical experiments (Sections 12.2.41 
andlH. 



1.4 Notation 



We gather here some useful notation and results. Let E z be the submanifold of W 
of positions at a fixed value of the reaction coordinate: 

E : = {qeR"^(q)=z}. (7) 

Let us introduce which is the probability measure fi conditioned at a fixed value 
of the reaction coordinate: 

expH^IV^-'^ 

"MX. = -J , (°) 

y^expc-js^iv^- 1 ^ 

where the measure Ox_ is the Lebesgue measure on E z induced by the Lebesgue 
measure in the ambient Euclidean space R" Z> E z . By construction, if Q is distributed 
according to the Gibbs measure ([TJ, then the law of Q conditioned to a fixed value 
z of ^ (Q) is iiz-- The measure | \~ l dOz, is sometimes denoted by 8^ q ^_ z {dq) in 
the literature. 

We recall the following expressions for the free energy A and its derivative A', 
also called the mean force (see 0): 



A(z) = -P -Mn^ Z-'expC-^IV^-^atJ , (9) 
A'(z) = f Fduz,, ( 10 ) 



where F is the so-called ZocaZ mean force: 



|v?| 2 r Viv? 

In the particular case when the reaction coordinate is just one of the cartesian coor- 
dinate, say 4 (g) = q", then 

A(z) = -p- l \n( [ Z- 1 exp(-pV(q 1 ,...,q"-\z))dq 1 ...dq"- 1 
\Jr"- [ 

and the local mean force is just F = d q "V, so that 

J M n-id q nV(q 1 ,...,q"-\z)ex V (-pV(q 1 ,...,q"- 1 ,z))dq 1 ...dq"- 1 



A'(z) 



Jk-i expC-jSV^ 1 , . . . )9 "-i,z))£V . ..</«?"- 1 
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2 Computing macroscopic stress-strain relations for 
one-dimensional chains of atoms 

In this section, we wish to compute the stress-strain relation of a one-dimensional 
chain of atoms, in the thermodynamic limit. More precisely, we consider a chain of 
I + N atoms, with its left-end atom fixed, and either submit the right-end atom to a 
force, and compute the average elongation, or prescribe the elongation, and compute 
the force. We will show that, in the limit N — s- °°, these two relations are identical, 
and that they can be computed in an extremely efficient manner. In short, passing to 
the limit N — >• °° makes tractable a computation that is, for finite and large N, very 
expensive. 

The relation between that question and the question of determining the free en- 
ergy of the system, when the reaction coordinate is the length of the system, will 
also be discussed. 

In the sequel, we first proceed with the nearest neighbour case (see Section lXTT ). 
We next address the next-to-nearest neighbour case in Section [Z2l which is techni- 
cally more involved. 



2.1 The nearest neighbour (NN) case 

We consider a one-dimensional chain of atoms, with positions q°, q l , . . . , q N . In this 
section, we only consider nearest neighbour interaction. In addition to this internal 
interaction, we assume that the atom at the right boundary of the chain is submitted 
to an external force /, and that the atom at the left boundary is fixed: q° = 0. The 
energy of the chain thus reads 



In the sequel, we will consider the limit when the number of atoms goes to °o. We 
wish to make sure that, even when N — > °°, the system occupies, on average, a finite 
length. To this aim, we introduce the rescaled positions u' = hq', with h = 1 /N. The 
energy now reads 



where again u = 0. 

For any observable <J>, depending on the variables u l ,...,u N , we define the 
canonical average of <P by 



2V 



E f (q 1 ,...,q N )=^W{q i -q i - 1 )-fq' 



i=\ 




(12) 




(13) 
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where the partition function Z reads 




exp (— fiEf (u l ,...,U 



A' 



)) du l ... du 



We assume in the sequel that W(r) grows fast enough to °° when |r| — » oo, so that Z 
is well defined (it is for instance enough that W(r) ~| r |^oo \r\ a with a > 1). 

We will be interested in the limit of ( < t')f i , when /V — > °°, and when <J> only 
depends on u N : <P(u l ,. . . ,u N ) = A(u N ) for a given function A. 

Remark 1. In ( fT3l l. we let the variables u' vary on the whole real line. We do not 
constrain them to obey u'~ l < u', which would encode the fact that nearest neigh- 
bours remain nearest neighbours. The argument provided here carries through when 
this constraint is accounted for: we just need to replace the interaction potential W 



2.1.1 Computing the strain for a given stress 

We first show a simple adaptation of [3] Theorem 1], which is useful to compute 
averages of general observables, in the thermodynamic limit, for the canonical en- 
semble at a fixed stress: 

Lemma 1. Assume that A continuous, that for some p>l, there exists a 

constant C such that 



by 




W(y) when y > 
+°° otherwise. 



Vy€R, \A(y)\<C(l + \y\P), 



and that 




Then 



\im (A(u N )) f N =A(y*(f)), 



with 



f R yexp(~p{W(y)-fy])dy 
k^P(-p[W(y)-fy})dy ■ 



(14) 



Proof. We observe that 
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where Wf(x) = W(x) — fx. Introducing y = , a change of variables in the 

h 



above integral yields 



« = z ~ l l A (jf i/j exp (-p f> w) ^ ■ ■ ■ dyN 

where, with a slight abuse of notation, Z = / exp I — /3 Wf (/) I ^y 1 . . . dy 



1=1 



Consider now a sequence |y } =1 of independent random variables, sharing the 
same law z~ l exp (—fiWf(y)) dy with z = / exp (—j3W/-(y)) c/y. It is clear that 



(A) f N = E 



i w ■ 

The law of large numbers readily yields that — ^ Y' converges almost surely to 

y*(f) defined by (0. 

We infer from [3] Theorem 1] that, for any force /, and for any observable A 
sufficiently smooth, the limit when — > °° of (A)j^ is 

Mm{A)^=A(y*(/)). 

Rates of convergence are also provided in the same theorem. □ 

Numerical simulations illustrating this result are reported in J3] Section 2.3]. 

In the specific case of interest here, namely computing the stress-strain relation, 
we take A(u N ) = u N , thus £#(/) := (A)C represents the average length of the chain, 
for a prescribed force /. We infer from the previous result that 

lim %(/) =y*(f). 

Af->oo 

We hence have determined the macroscopic elongation, namely y*(f), for a pre- 
scribed microscopic force / in the chain. 

Notice that, in this specific case, A is a linear function, so we actually have 
£jv(/) = y*(f) for any N. The result of Lemma[T]remains interesting for computing 
standard deviation of the average length, for example. 

Remark 2. The force between atoms j and / — 1 is W . Its canonical 

V h J 

average, defined by ( fl3l l, is 

C ] N = Z~ X ( W' f uJ ~ uJ \exp(-pE f (u l ,...,u N ))du 1 ...du N 
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= z 1 W (y ) exp (-pjt[W (/') - fy'] jdy l ...dy N 

Jg W (y ) exp (-Ji [W (y ) - fyi] ) dy> 
J R exp(-p{W(yj)~fyJ})dyj 

Jg [w (y ) - /] exp (-p [w (y ) - /y] ) </y 

f M exv(-P[W(yJ)-fyj])dyJ 

, u j -u j - 1 

where > J = . Integrating by parts, we see that the second term of the 

h 

last line vanishes. We hence obtain that the average force between two consecu- 
tive atoms is independent of j (the stress is homogeneous in the material), and is 
equal to its prescribed microscopic value /: 

v,\ v#, 4=f. 

Imposing a force / on the right boundary atom hence implies that the average force 
between any two consecutive atoms is equal to /. 



2.1.2 Computing the stress for a given strain 

In the previous section, we have prescribed a force, and computed an average elon- 
gation. We now prescribe the length of the material, by imposing «° = and u N = x 
(see Fig . |TJ . 



N 

Fig. 1 One-dimensional chain of 1 + N atoms, where the total length of the system is prescribed 
at the value x. 



As we fix the position of atom N, the system is insensitive to any force / imposed 
on that atom. We hence set / = 0. Our aim is to compute the force in the chain, 

f W'( X ~ uN l ) exp (-pE (u l , . . . , u N -\x)) du l ...du N - 1 

W= kN ~ l } k J , (15) 

/ exp (— pEo(u l , . . . ,u N ,x)) du x ...du N 1 



or, more precisely, its limit when — >• °°. Note that, as all the play the 

h 

same role in the above expression, we also have, for any 1 < i < N — 1, 
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i-1 ■ 



W' [ ^—^ — J exp(-j8£ (" 1 ,---,M A ' _1 ^)) du 1 ...du 
exp(-j8£ (« 1 ,---,M Ar_1 ^)) du 1 ...du N - x 



sr N (x) 



The force between atom N and AT — 1 is thus equal to the force between any two 
consecutive atoms. 

We infer from (fT~5T > that ,%i(x) = F^(x), where 

Fn(x) = — — — In f exp (— f}Eo(u l , . . . 7 u N ~ l ,xj) du l ...du N ~ l 
pN [Jm n - 1 

Hence NFn is the free energy of the material associated to the reaction coordinate 
% (u , . . . , it) = u , and Fn is a rescaled free energy (free energy per integrated out 

u' — u'~ 1 

particle). Using the variables y' = , we also see that exp(— fiNFtj{x))dx 

h 

is (up to a normalizing multiplicative constant) the probability distribution of the 

1 N N 

random variable — V Y\ when is a sequence of independent random vari- 

ables, sharing the same law exp (~fiW{y))dy, with z = / exp(— [3W(y))dy. 

Jr 

1 , 

In the case W(y) = —(y — a) , it is possible to analytically compute Fn(x), and 

to observe that there exists a constant Cjy, independent of x, such that Fn(x) + Cn 
has a finite limit when N — > °°. In the general case, the limit of Fn is given by the 
following result, which relies on a large deviations result for i.i.d. random variables: 

Lemma 2 ([31, Theorem 2). Assume that the potential W satisfies 
V^eK, [ ex.p(£y-pw(y))dy < +°°, 



and exp(-j3W) e H l (R). Then 

J3 AT 



lim ( F N {x) + - In - 1 = F„(x) (16) 



with 

F oa (x) := — sup ( %x — In 

P 56 



z- 1 / exp(^-j3W(y))^ 



(17) 



ana? z= exp(—f}W(y))dy. This convergence holds pointwise in x, and also in 

Jm 

ZZ , for any 1 < p < °o. As a consequence, F N converges to F^ in ' p . 

We hence obtain the macroscopic force Fi,(x) for a prescribed elongation x. Nu- 
merical simulations that illustrate this result are reported in [3, Section 2.3]. 
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Remark 3. The additive term /J -1 \n(z/N) in ( TT~6b can be seen as a normalizing con- 
stant. Indeed, as mentioned above, NFn is a free energy, and the correct normaliza- 
tion for exp(— fiNFfj) to be a probability density function is: 

dx= 1. 



Remark 4. Fn is a challenging quantity to compute. One possible method is to com- 
pute, for each x, its derivative FL(x), and deduce fjy (this is the so-called thermo- 
dynamic integration method). Note that F^(x) = !J^(x) is given by ( fT3T >: it is a 
canonical average of some observable, in a space of dimension N — 1 ^s> 1 . In con- 
trast, Foo is easier to compute, since it only involves one-dimensional integrals or 
optimization problems, o 



/ 



F N (x) + jln^ 



2.1.3 Equivalence of stress-strain relations in the thermodynamic limit 

The function we maximize in ( fTTI i is concave, so there exists a unique maximizer 
4 (x) in (fTTT i. that satisfies the Euler-Lagrange equation 

_ f R yexp(Z(x)y-PW(y))dy 
- ! R e W ($(x)y-pw(y))dy ' ' 



We observe that 



On the other hand, recall the definition (fT~4T > of y*(f): 

Lycxp(-p[W(y)-fy])dy 



f{f) 



/ K exp(-/3 [W(y)-fy})dy 



Comparing ( fT~8T > and ( fPfl i. we see that y* (j3 1 <^ (x)) =y*(F^(x)) — x. The function 
/ 1 — ^ y*(f) is increasing (because its derivative is positive), thus it is injective, and 
we also get the converse relation: F^,(y*(f)) = f. 

Otherwise stated, the relation / h-> y*(f) and x H> Fl,(x) are inverse one to each 
other. So, prescribing a microscopic force / and computing the macroscopic elon- 
gation is equivalent to prescribing an elongation and computing the macroscopic 
force, in the thermodynamic limit (namely in the limit N — > °°). 



2.2 The next-to-nearest neighbour (NNN) case 



We now consider next-to-nearest neighbour interactions in the chain. Again, the first 
atom is fixed: u° = 0, whereas the last one is submitted to an external force /. The 
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(rescaled) energy reads 



*>v ■") \r-r- ) + g * j -/i- (19) 

If W2 = 0, this energy reduces to (1121 1. Averages of observables are again defined 
by ([T3ll. 



2.2.1 Computing the strain for a given stress 

Our aim, as in Section 12.1.11 is to compute the macroscopic strain, which is the 
average length of the material, that is 



%(/) = <"% 



where (-) N is the average with respect to the canonical measure associated to Er. 
We introduce the notation 

W lf (x)=W l (x)-fx, 

which will be useful in the sequel. A simple adaptation of |f3] Theorem 3] yields the 
following general result: 

Lemma 3. Assume that A : R 1— > R is continuous, and that there exists p > 1 and 
C > such that 

\A(x)\<C(l + \x\P). 

Assume also that Wif and W2 both belong to L\ c (R), that they are bounded from 
below, and that, for any x £ R, we have |Wi/(jc)| < 00 and \W2(x)\ < °°. In addition, 
we assume that e~P w 'f ande~P Wl both belong to (R), with 

[ (1 + \x\") e'Wiftodx < +00 and f (1 + \x\ p ) e~ m{ ^dx < +<*>. 

JR JR 

\im(A(u N )) f N =A(y*(f)) (20) 



Then 

N 



with 

/(/)= / yvj(y)dy, (2i) 

where Xfff solves the variational problem 

X f = max ( f y,(y) V(z) K f (y,z) dydz; j y 2 {y)dy = l) , (22) 
yeL 2 (R) [Jr 1 Jr J 

with 
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K f (x,y) 



exp 



-pW 2 (x+y) - |wi f ( x ) - |wi ,(y) 



(23) 



We only provide here the main arguments to prove this result (see [3] Sec. 3.1.1 
and Theorem 3] for details). They will be useful in the sequel. The observable A (u N ) 
only depends on ur , thus 

(A(u N )) f N = Z 1 [ A(u N )exp(-l3E f (u\...,u N ))du 1 ...du N 



N-l 



u i+1 -i/- v 



du ... du N . 



Introducing again the variables y' 

(A(u N )) f N =Z- 1 



-, we see that 



J MN A U tA «p i-PWif (y 1 )) ft*/ (^V) ■■ ■ d f 

(24) 



with 



k f (y-\y') = exp (-pWif (y) - j3W 2 (y^ 1 +/')) 
Assume for a moment that / kf{a,b)db = 1 . Then we see that 



(A(« A, )>^ 



where {F'|. =1 is a realization of a Markov chain of transition kernel kf, and where 
T 1 has the initial law (up to a normalization constant) exp(— fiWif (y 1 )) dy 1 . A 
law of large numbers argument, now for Markov chains, yields the large limit of 
(A(u N ))-jy (recall that, in the case of the NN model considered in Section l2.1.1l this 
limit is given by a law of large numbers argument for i.i.d. sequences). 

In general, of course, / kf(a,b)db ^ 1 . There is thus a slight technical difficulty 

m 

in identifying a Markov chain structure in (I24l i. It yet turns out that the above argu- 
ment can be made rigorous as follows. Consider the variational problem d22l , with 
Kf defined by (1231 . Under our assumptions, Kf £L 2 (RxR). Using standard tools 
of spectral theory of self-adjoint operators (see e.g. ifTOl ). one can prove that this 
problem has a maximizer (denoted i//y), and that, up to changing \\ff in — V/, the 
maximizer is unique. In addition, one can choose it such that \\ff > 0. We can next 
define 



(25) 
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which satisfies 

/ g f (y,z)dz = l, f V}(y) 8f(y,z)dy = \l/j(z). 
Jr Jr 

The average d24b now reads 



{A{^%=Z? L^{jff/j Wfiy')e-^f^ 



*gf(y l S) • ..gf(y N -\y N ) ^ V ■ ■ ■ *? 



(26) 



withZ, = Wf (y l ) e "!*W) gf (y\yl) . ..g f (f-\f) e 
Thus 



where (T 1 , . . . , T") may now be seen as a realization of a normalized Markov chain 
of kernel gf, with invariant probability measure \\fj. 

Under our assumptions, the Markov chain has a unique invariant measure, and 
satisfies a law of large numbers with respect to it. This yields the convergence ( f20b . 
Numerical simulations illustrating this result are reported in J3J Section 3.1.3]. 

In the specific case of interest here, namely computing the stress-strain relation, 
we take A(u N ) = u N , thus £#(/) := (A)C represents the average length of the chain, 
for a prescribed force /. We infer from the previous result that 

lim e N (f) =y*(f). 

We hence have determined the macroscopic elongation, namely y*(f), for a pre- 
scribed microscopic force / in the chain. 

We conclude this section by showing the following result, which will be useful 
in the sequel. 

Lemma 4. Under the assumptions of Lemma \3\ introduce the asymptotic vari- 
ance a 1 (/) defined by 

o\f) = [ (x-y*(f)) 2 Y#*)«fr + 2£E((?(-y*Cf))(?i -/(/))) (27) 

JR ; > 2 V ' 

where (y^J is a Markov chain of transition kernel gf, and of initial law \\Pj, the 
invariant measure. 

Assume that 2 {f) ^ almost everywhere. Then the function f i— > y*(f) is in- 
creasing. 
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Note that the right-hand side of (|27T i is exactly the variance appearing in the Central 
Limit Theorem for Markov chains 11221 Theorem 17.0.1]. It is thus non-negative. 

More precisely, we have that lim NYar Y^j = 2 (f) where (Yij is the 

Markov chain defined in the above lemma. 

Proof. Let £#(/) := (u N )^. An analytical computation shows that 

(("" ) 



,N\2\.f f/„N\f^ : 



' N 



Thus the function / n- £#(/) is non-decreasing. By Lemma[3] y*(f) is the pointwise 
limit of £/v(/) : ^ i s mus non-decreasing. It remains to prove that it is increasing. 

Let us now compute the limit when N —> °° of Djy(/). Using Theorem 4], we 
see that 

lim D N (f) =J3 a 2 (/), 

2V->oo 

where <J 2 (f) is defined by d27l i. 

Let us now fix T and T > T. Since Dn(/) > 0, we can use Fatou lemma, which 
yields that 

J3 r a 2 (f)df= f\imMD N (f)df<limmf [* D N {f)df =y*(T)-y*(z). 

Jz Jz Jz 

As <7 2 (/) > almost everywhere, we thus obtain that x i-» is an increasing 
function. □ 



2.2.2 Computing the stress for a given strain 

We now prescribe the length of the material, by imposing u° = and u N = x. Our 
aim is to compute the average force in the chain, 

/ A h (u N -\u N - 2 -x)^v(-pE Q (u\...,u N -\x)) du l ...du N - 1 

y N t x \ = Mil ; 

[ exp(-BE (u\...,u N - 1 ,x)) du l ...du N - 1 

(28) 

where Eq is the energy (fT9b with / = 0, and where the observable A/, is the force 
acting at the end of the chain, which reads 

A h (u N -\u N - 2 ;x) = W{ (^^) +Wi (^^) • 

More precisely, we are interested in lim ,%i(x). 

As in Sec tion l2. 1.21 we see that &n(x) = F^(x), with 
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exT>(-BE (u 1 ,...,u N - 1 ,x)) du\..du N - } 



(29) 



Again, NF^ is the free energy associated to the reaction coordinate t, (u l , . . . , u N ) = 
ur , and Fn is a rescaled free energy (free energy per integrated out particle). In 
the NN case, we have computed the large N limit of Fn(x) using a large deviations 
result for i.i.d. random variables. Comparing Sections |2. 1 . 1 1 and 12.2.11 we also see 
that moving from a NN setting to a NNN setting implies moving from a framework 
where random variables are i.i.d. to a framework where they are a realization of a 
Markov chain. It is hence natural to try and use a large deviations result for Markov 
chains to compute the large N limit of (|29l . 

We now assume that the underlying Markov chain satisfies the following point- 
wise large deviations result: 

Assumption 1 Consider the Markov chain {J" } of kernel k £ L 2 (lx R). As- 
sume that, for any B, £ R, the function exp(^y)k(x 7 y) £ L 2 (R x R). 
Introduce the operator (on L? 1 



(Qt < P)(y)= / (p(x)exp(£,y)k(x,y)dx 



and assume that it has a simple and isolated largest eigenvalue A ), and that 
4 i— > In A ) is convex. 

- 1 N ■ 

Let exp(— NFn{x))dx be the law of the random variable — J^'. We assume the 



large deviations principle 



N:._ A 



lim F N (x) =Foo(x) (30) 



where 

Foo(jc):=sup(§JC-lnA(§)). (31) 

We moreover assume that the convergence holds pointwise in x, and also in 
Lf , for any 1 < p < °o. As a consequence, f' n converges to in W lo ^' p . 

Note that similar results in a finite state Markov chain setting are reviewed in 15] 
pages 60-61] or (8] Sec. 3.1.1] (the continuous state case is addressed in e.g. (8] 
Sees. 6.3 and 6.5]). In the discrete state case, one can prove that ^ i-> lnA(^) is 
convex (see j9j Exercise V.14]). We will numerically check in the sequel that this 
assumption is indeed satisfied in the example we consider (see Fig.|2]i. 

Remark 5. We have assumed that the operator has a simple and isolated largest 
eigenvalue. This can be proved for many kernels k, using for instance Krein-Rutman 
theorem 11251 . In the case of interest in this contribution, we will use the specific 
expression of the kernel to transform the operator into a self-adjoint Hilbert- 
Schmidt operator on L 2 (R) (see Remark [7] below). We will thus be in position to 
work with self-adjoint compact operators, o 
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Remark 6. In the NN case, when k(x,y) = 6(y) = z exp(— j3W(y)), the sequence 
{F'I^j is a sequence of i.i.d. variables sharing the same law 6(y)dy. The operator 

Qf has a unique eigenvalue A (B) = / exp(^y) 9(y)dy. We then recover the large 
deviations result of i.i.d. sequence given in Lemma[2](see e.g. llT^[T3lfT4ll26l '). o 

We now wish to use Assumption Q] to compute the large N limit of (19) . As 
pointed out in Section 12.2.11 there is a slight technical difficulty in identifying a 
Markov chain structure in the NNN setting, related to the normalization of the 
Markov chain kernel. We thus cannot readily use Assumption Q] We now detail 
how to overcome this difficulty. 

Consider an observable A that depends only on un. In view of (|29| | and (|26| |. its 
canonical average reads 

(AW=Z- 1 / A (u N ) exp (~f$Eo (u 1 ,.. . ,u N ~ l ,u N )) du l . ..du N ~ l du N 

= Z _1 [ A(x)exp(-PNF N (x))dx 



XgoW J Z ) ■ • -goif-'y) ^Tyv^ dy l . ..dy N , 

where go is defined by d25l l and x/zq is the maximizer in (122b . when the body force / = 
0. Let ^(y 1 , ■ ■ ■ ,y N ) be the probability density of a Markov chain {l"}.^ of ker- 
nel go, where the law of Y 1 is (up to a normalization constant) XjfQ (y l ) e~ % Wl & ' dy l . 
Then 

/ A(x)exp(-pNF N (x))dx = C N [ a( if y' ) &>(y\ ...,y N ) r(y N )dy l .. .df 

«/ffi JWL \ • -J J 

(32) 

where Cn is a constant that does not depend on the observable A, and 

rif) 




Wo(y N ) ' 

Let now a N (x,y N )dxdy N be the law of the couple I j^J^Y l ,Y N 1 . We recast $32% 
as 

A(x) exp (-PNF N (x)) dx = C N [ A (x) a N (x,y N ) r (y N ) dxdf . 

Jm. 2 



As this relation holds for any observable A, with a constant Cn independent of A, 
we obtain 

cxp(-pNF N (x)) =C N [ a N (x,y N ) r{f) dy N . 
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Assuming that r and 1 jr are in L°°(M), we have 

Qvll l/r\\ L l [ a N (x,y N ) dy N < exp (-pNF N (x)) < C N \\r\\ L - [ a N ( Xl y N ) dy N . 
Jr Jr 

As a consequence, since the function r is independent of N, 

k ln i aN(xy)dy 



lim (Fn(x) +Dn) = lim 

N— ><*• N— ><»= 



(33) 



where Dm = — — In Civ. Recall now that 
/3iV 



Yn(x) = a N (x,y N ) dy N 
Jr 

1 N N 

the density of — where {F'} ( ._j is a realization of the Markov chain of 



l 

kernel go- The behaviour of Jn when N — > °° is given by Assumption[T] 

Urn -~\nyb(x)=T„{x), (34) 

jV->+<*> N 

where Foo is given by (T3TT >. Collecting d33l l and ( f34t . we hence obtain that 

lim (F N (x)+D N ) = -F„{x). 
We thus have the following result: 

Lemma 5. Assume that W\ and W2 both belong to Lj (R), that they are bounded 
from below, and that, for any x £ R, we have \Wi(x)\ < °° and \W2(x)\ < °°. In 
addition, we assume that e-P w > and e~P W2 both belong to W^ 1 (R), with 

f e -Pw 1 (x) dx<+00 and f e -W4 <+0O] 
Jr Jr 

and that, for any £, £ R, we have exp {E,x - pW\ (x)) £ L l (R). 

Under Assumption \l\for the kernel go defined by ( 1251 ), the limit of \29\ is given 

by 

lim {F N (x) + C N )=F 00 {x) (35) 

where Cn is a constant that does not depend on x, and Fee is given by the Legendre 
transform 

F»(*):=-isup(£c-lnA(§)) (36) 

P $eR 

where A(^) is the largest eigenvalue of the operator ( defined on U 
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(Qt<P)(y)= / 9(x)exp(&)go(x,y)dx (37) 

where go is defined by ( 125P . 77ie convergence ( I55P /io/ofs pointwise in x, and also 
in , for any 1 < p < °o. As a consequence, the macroscopic force in the chain 

&n(x) = Flf(x) converges to F^ in W loc ' p . 

We hence obtain the macroscopic force F^ (x) for a prescribed elongation x. Note 
that, under our assumptions, in view of its definition d36b . F^ is (up to the factor 
/3) the Legendre transform of some function. It is hence always a convex function. 
Thus, as in the zero temperature case, we observe, in this one-dimensional setting, 
that the macroscopic constitutive law x n- is a convex function. 



Remark 7. In view of the definition (1251 1 of go, we see that 

— ^—r\ — = 7- / — tt sx P(hy)Ko(x,y)dx. 
Thus A (4 ) is also the largest eigenvalue of the operator 

(Q$P)Cy) = V" / ?(*) «p(§y)Jeb(jc,y)<fa. 
Furthermore, if A is an eigenvalue of , then 

<p(x) exp(%y)K (x,y)dx = A A<p(y) 



where <p is an associated eigenfunction. Thus 

exp^y/2)exp(^x/2)Ko(x,y)dx = XoX- 



exp(^/2) -r^,-,-™-,-,"^- «p(§y/2) 
and AoA is an eigenvalue of the operator 



(G£9)Cy)= / <pWexp(^y/2)exp(^x/2)^o(^y)^- 

The converse is also true. As A(^) is the largest eigenvalue of the operator Q^, we 
have that AqA {£, ) is the largest eigenvalue of the operator . 

Note that is a self-adjoint compact operator on L 2 (R), which is thus easier 
to manipulate theoretically and numerically than . In particular, using standard 
tools of spectral theory of self-adjoint operators (see e.g. iflOl ). one can prove that 
the largest eigenvalue of is simple, and that the associated eigenvector (which 
is unique up to a multiplicative constant) can be chosen such that > 0. o 
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2.2.3 Equivalence of stress-strain relations in the thermodynamic limit 

In Section 12.2.11 we have identified the function / H> y*(f), that associates to a 
prescribed force / the macroscopic elongation y*(f). Next, in Section 12.2.21 we 
have identified the function x i— > Fl,(x), that associates to a prescribed elongation x 
the macroscopic force F^(x). We show now that these functions are reciprocal one 
to each other. 

Consider the optimization problem d36l l. Since the function £ n- In A (£, ) is con- 
vex (see Assumption[T|i, there exists a unique maximizer % (x) in (l36l l. which satisfies 
the Euler-Lagrange equation 



We also observe that 




We see from (l38l l that we need to compute A'(§). Recall that A (£, ) is the largest 
eigenvalue of the operator (1371 . In view of Remark [7] AoA(^) is also the largest 
eigenvalue of Q^. Denoting the associated eigenfunction satisfying ||f^|| L 2 = 1 
and > 0, we thus have 




where 



K${t,y) = exp(<^/2)exp(^/2)K (f,y). 
Multiplying by 5^ (y) and integrating, we obtain 



(39) 




(40) 



We thus have, using that is self-adjoint, 







(41) 
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Collecting (|38j, (|40]l and (ED, we see that 



,^ w Cy)^w(04 w (*oO*<'y 
y9% w Cy)^w(0«o W ('.3')*«' 



,^ w (y)^ w (04 w ('.y)* d y 



(42) 



where we have used, at the second line, flu*Jsf«M 

On the other hand, we have obtained that the macroscopic elongation y*(f), for 
a prescribed force /, is given by (f2TT >. namely 



y*(f)= I yv}(y)dy 



where is the maximizer of the variational problem d22l . As Kf is symmetric, the 
Euler-Lagrange equation of $2% reads 



I 

J 7: 



\lf f (t)Ko(t,y)ex V \Bf?±y) dt 
Vf (t)K^ f (t,y)dt 

where is defined by (f39l >. Thus y/f is an eigenfunction associated to the largest 
eigenvalue A/ of the Hilbert-Schmidt operator Qpf of kernel Kq . By definition of 
Wpf, and using the fact that the largest eigenvalue of Qp / is simple, we obtain 

k 

V 



¥p f = ±Yf and AQ3/) = ^. 



Thus 

f{f)= I yVp f (y)dy. (43) 
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We deduce from the comparison of (@2]i and (|43j that y* (j3 "~ 1 <fj (x) ) = (i^ (jc) ) = x 
Recall now that the function / H> y*(f) is increasing, as shown by Lemma|4] It is 
thus injective, and we also get the converse relation Fl,(y*(f)) = f. 

As a consequence, as in the NN setting considered in Section l2.1.3l the relation 
/ ^ y*if) an d x ^ FL{x) are inverse one to each other. Prescribing a microscopic 
force / and computing the macroscopic elongation is equivalent to prescribing an 
elongation and computing the macroscopic force, in the thermodynamic limit. 

2.2.4 Numerical computation of F^ and comparison with the zero 
temperature model 

For our numerical tests, we follow the choices made in j'3||, for the sake of compari- 
son. We thus take the pair interaction potentials 

Wl (x) = ^(x-l) 4 +±x 2 and W 2 {x)= l -{x~2.\)\ 

Note that these potentials satisfy all the assumptions that we have made above. 

We are going to compare the free energy derivative ,%i(x) = F^{x) with its ther- 
modynamic limit approximation F^(x). The reference value F^(x) is computed as 
the ensemble average (|28l l, along the lines of Q-©. To compute we proceed 

as follows: 

(i) We first compute the largest eigenvalue A(^) of the operator 1371 . for all B, in 
some prescribed interval. 

(ii) For any fixed x in a prescribed interval, we next consider the variational prob- 

lem (l36l l. compute its maximizer £, (x), and obtain F^,(x) using F^(x) = . 

r 

In practice, using Remark[7] we work with the operator Q^, which is easier to ma- 
nipulate since it is self-adjoint and we do not need to first solve (f22b . We thus first 
compute the largest eigenvalue XqA (£, ) of , and next compute the Legendre trans- 
form of the function £, n- ln(AoA (£,)). The maximizer is the same as that for foo(x). 
On Fig. |2] we plot the function £, M> In (AoA (<!;)), and observe that it is convex, in 
agreement with AssumptionQ] 

We first study the convergence of F^(x) to Fl,(x) as N increases, for a fixed chain 
length x = 1 .4 and a fixed temperature 1/J3 = 1. Results are shown on Figure[3] We 
indeed observe that F^(x) — > F^(x) when N — ^ +°°. 

We now compare F^(x) with its approximation F^ix), for N = 100 and = 1. 
Results are shown on Figure|4] We observe that F^(x) is a very good approximation 
of Ffl{x), for any x in the considered interval. 

For the sake of comparison, we now identify the zero temperature behaviour of 
the system, in the thermodynamic limit. At zero temperature, for a finite N, we 
model the system by minimizing the energy £0, with prescribed Dirichlet boundary 
conditions (this corresponds to prescribing the elongation, and computing the force; 
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Fig. 2 Plot of ln(AoA(£)) as a function of % (temperature 1/j8 = 1). 
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Fig. 3 Convergence of (x) (shown with error bars computed from 40 independent realizations) 
to f^(jc) as N increases (temperature 1//3 = 1, fixed chain length x = 1.4). 



alternatively, one could impose Neumann boundary conditions, i.e. prescribe a force 
and compute an elongation): 

J N (x) = — inf {E Q (u°,u\ . . . , u N - 1 , u N ) , u° = 0, u N = x) . (44) 

We have the following result, which proof will be given below: 
Lemma 6. Let us introduce <p defined by 

<j)(x)=Wi(x)+W 2 (2x). (45) 

Assume that there exists a > such that 

Wi(x) > ax 2 , (46) 
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and that W\ and are non-negative and strictly convex functions. Then we have the 
pointwise convergence 

lim Jn(x) =<j)(x). 

Assume in addition that <p £ L^ oc for some 1 < p < °° and that W2 is non-negative. 
Then the above convergence also holds in £| oc - As ci consequence, J^(pcj converges 
to <j)'(x) inW^f. 

When the temperature is set to zero, the energy thus converges, in the thermo- 
dynamic limit, to (j)(x), and the force (i.e. the derivative of the energy with respect 
to the prescribed Dirichlet boundary condition) converges to (j)'(x). We plot on Fig- 
ure|4]the function* 1— > <j>'(x). We clearly observe the effect of temperature, as F^(x) 
for j3 = 1 significantly differs from 0'(x). 




Fig. 4 We plot F^(x) and Fi(x) for the temperature 1/j8 = 1 and N = 100. On the scale of the 
figure, Ffl(x) and F^(x) are on top of each other. We also plot the zero temperature response <j)'(x). 



Proof (Lemma\6ty. Let 

X N (x) = {(u°,u 1 ,...,u N - 1 ,u N )eR 1+N , m° = 0, u N =x} 

be the variational ensemble for the problem ( l44b . The configuration u' = ix/N 
clearly belongs to that ensemble. We thus obtain the upper-bound 

Jn(x) < Wi (x) + ^-W 2 (2x). (47) 

In the sequel, we first show a lower-bound for Jn(x), and next study its behaviour 
when N — > °°. 

Let us first build a lower bound for Jn(x). Assuming for the sake of simplicity 

that is even, and using the short-hand notation y = , we have 

h 
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N 



iv (=1 

= ^Cy 1 ) + ^i(y v ) + ^ f[wi (f'-^+w, (/)] 



2JV 



2JV 



, N/2-l 

\f E [wiC^+^C^ 1 )] 



2N 



i=i 



By convexity of W\, we obtain 



1 w 



u' — u l 1 



1 1 l"/ 2 



2JV 



N 



N/2-1 



N ' 



Taking into account the next-to-nearest interactions, we thus obtain that, for any 

( a °y,..,/-7)ei M , 



2 ,V 



,.2i-l , ,.2<\ 



N/2-l 



+- E 



i=\ 



(y 2i +y 2i+1 ) 



where is defined by d45t . As is convex, we deduce that 



I 1 1 1 / 1 N/2 

-Eo(u) > — W l (y l ) + — Wi(v N ) + -(j) -E W 
N 2N yJ ' 2N v 1 2 Y \ N u L " 



i=\ 



N-2 
2N 



( 1 N > 2 - 1 \ 



^i(y)+^w 1 (/)+^(^- M °) 



N-2 



'-<$> 



N 



2N \N-2 

As a consequence, for any configuration u G X^{x), we have 



1 



Eof" ," 1 ,...,^ 1 ,^) >'E N (u i y- l ',x) 



(48) 



with 



E N (u l ,J-\x) = ^W x {Nu l ) + ^W x {N{x-u N - x ))+ l -<$>{x) 
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, N-2 ( N , N _i x , 

H <M ( u —u) 

IN Y \N-2 K > 

We infer from d48b the lower bound 

J N (x) > J N (x) (49) 

with 

7 N (x) =mt{E N (u\u N ~ l ;x); u l £ R, u N ~ l e R} . (50) 

We now study the auxiliary variational problem (l50l to determine the limit of 
Jn(x) when N —> °°. Since (j) is non-negative, we infer from (|46| | that 

E N {u\u N - l -x) > ^ [(u l ) 2 + (x-u N - 1 ) 2 } > 0. 

As a consequence, Jn(x) > 0, and any minimizing sequence is bounded. Up to ex- 
traction, it thus converges to a minimizer, that we denote (m 1 ,^ -1 ). As W\ and 
are strictly convex, it is easy to see that the hessian matrix of En is positive definite, 
hence En is also strictly convex, hence it has a unique minimizer. The problem ( T50b 
is thus well-posed. To underline the dependency of its minimizer with N, we denote 
it (m 1 (N), u n - x (AO) in the sequel. 

The Euler-Lagrange equation associated to d50b reads 

W[{Nu\N)) = f (J^ (u N - l (N)~u\N))^j =W{(N(x-u N - 1 (N))). 

As W\ is strictly convex, this implies that 

u l (N) =x-u N - 1 (N), 

N , _ (51) 



where the function % = (W^) -1 o 0' is independent of N, and increasing. 

Let us now show that u 1 (N) is bounded with respect to N. If this is not the case, 
then, without loss of generality, it is possible to find a subsequence (p(N) such that 
limyv^oo m 1 (q>(N)) = +°°. Passing to the limit in the second line of (fSTb . one obtains 
a contradiction. Thus u l (N) is bounded. 

In view of the first line of ( l5Tb . u N ~ l (N) is also bounded. Up to a subsequence 
extraction, (77 1 (N),u N ~ l (N)) converges when N — > °° to (u 1 (oo),u N (°°)). We infer 
from ( BTT ) that Ti l (oo) = Q and ^"'(qo) = x, thus the limit is unique, and all the 
sequence converges: 

lim u l (AO = 0, lim if' 1 (N) = x. (52) 



We next infer from the above limits and ( Bil l that 
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lim Nu l (N) = lim N(x-lf- 1 (N)) = Y(x). (53) 

By definition, we have 

J N (x) = mf{E N (u 1 ,u N ~ l \x); u l el, k^ -1 el} =E N (u l (N),u N ~ l (N);x). 
In view of d52l and ( f53l l. we obtain 

Mm J N (x) = lim E~n(u 1 (N),u n ~ 1 (N);x) =Q(x). (54) 

Collecting d47| i, d49l > and (|54l , we obtain the claimed pointwise convergence of 
Jn(x) to 0(x). 

Under the additional assumption that W2 is non-negative, we deduce from d47l i 
that, for any N and any x, 

< J N (x) < Wi (x) + W 2 (2x) = <j)(x). 

As € L[^ c , we obtain the convergence of Jn to <j) in L^ c . □ 



3 A coarse-graining procedure in the dynamical setting 

In this section, we present a procedure for coarse-graining a dynamics. More pre- 
cisely, we consider Q t G W solution to the overdamped dynamics (fjjl, and a reac- 
tion coordinate £ : R" n- R. Our aim is to find a closed one-dimensional dynamics 
of type (0 on a process rf t , such that rf t is a good approximation of £,(Qt). In 
Sections [3 .21 and [3~3l we build such a process (see (l65l l below), and present an an- 
alytical estimation of its accuracy (the obtained estimate is an upper-bound on the 
"distance" between the laws of £, (Q t ) and rf t at any time f). We will next report on 
some numerical experiments that somewhat check the accuracy of rf t in a stronger 
way (Section [3~4l . 



3.1 Measuring distances between probability measures 

We introduce here some tools that will be useful in the sequel, to measure how 
close two probability measures are. Consider two probability measures v(dq) and 
rj(dq). The distance between the two can be measured by the total variation norm 

|| v — TJ||tv. which amounts to the L'-norm / | V / v(?) — VriW) \ dq in case v and 77 
have respectively the densities y/ v and yfri with respect to the Lebesgue measure. 

When studying the long-time behaviour of solutions to PDEs (such as long time 
convergence of the solution of a Fokker-Planck equation to the stationary measure of 
the corresponding SDE), the notion of relative entropy turns out to be more useful. 
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Under the assumption that v is absolutely continuous with respect to 77 (denoted 
V 77 in the sequel), it is defined by 

H{ ?^=h{%) dv - 

The relative entropy provides an upper-bound on the total variation norm, by the 
Csiszar-Kullback inequality: 

||v-T7|| T v< V2#(v|f?). 

In the sequel, we will also use the Wasserstein distance with quadratic cost, which 
is defined, for any two probability measures v and 77 with support on a Riemannian 
manifold E, by 



W(v,rj) = J inf / d E (x,y) 2 n(dx,dy). (55) 
y 7ien(v,ri)JzxL 

In the above expression, d% (x,y) denotes the geodesic distance between x and y on 




d z (x,y) =inf { \l [ \a{t)\ 2 dt; aeC 1 ([0,l] ) I),a(0)=i,a(l)=j 



and n(v,Tj) denotes the set of coupling probability measures, that is probability 
measures n on E x E such that their marginals are v and rj: for any test function <J>, 

/ <P(x)n(dx,dy) = / <&(x)v(dx) and / &(y)n(dx,dy) = / <P{y)r\(dy). 
Jexe Je Jexe Je 

In the sequel, we will need two functional inequalities, that we now recall (TJ: 

Definition 1. A probability measure 77 satisfies a logarithmic Sobolev inequality 
with a constant p > if, for any probability measure v such that v 77, 

H(V|77) < ^/(V|T7) 

where the Fisher information /(v|?7) is defined by 

f (vi„)=/v b (|Y 2 



dv. 



Definition 2. A probability measure 77 satisfies a Talagrand inequality with a con- 
stant p > if, for any probability measure V, 



Some remarks on free energy and coarse-graining 



29 



W(v,J7)<y-//(v|77). 

We will also need the following important result (see 1231 Theorem 1] and (4)): 

Lemma 7. If rj satisfies a logarithmic Sobolev inequality with a constant p > 0, 
then r\ satisfies a Talagrand inequality with the same constant p > 0. 

The following standard result illustrates the usefulness of logarithmic Sobolev 
inequalities (we refer to [HE] [27] f° r more details on this subject). 

Theorem 1. Consider Q t solution to the overdamped Langevin equation (0, and 
assume the stationary measure y*>{q)dq = Z exp(— fiV(q))dq satisfies a loga- 
rithmic Sobolev inequality with a constant p > 0. Then the probability distribution 
\jf(t, •) of Q t converges to \jfoo exponentially fast, in the sense: 



Vf >0, //(y/(r,-)|y/oo) <//(v/(0,-)|v / »)exp(-2pj3" 1 f) 



(56) 



Conversely, if ( 1561 ) holds for any initial condition 1/(0, ■), then the stationary mea- 
sure V / o°(?) dq satisfies a logarithmic Sobolev inequality with a constant p > 0. 

Proof. The probability distribution function yf{t,q) of Q t satisfies the Fokker- 
Planck equation 

d,\j/ = drv OVVj+jfr'ziy/. (57) 
As Vy/oo = — j3 y/ooVV , we recast the above equation as 



d t \j/ = p l div 



y/oo 



Note that this equation implies that 
relative entropy 



<f(*)=#(y(v)lv~) = / i 



y/(t,q)dq is a constant. Introduce now the 

Y(f,qY 



\\r{t,q)dq. 



Then 



dS 
dt 



¥_ 
<//«, 

¥_ 



In 

! 

In 

! 

-1 

= -/3- 1 /( V /(f,-)l^ 



a f — - 1// 



J3 _1 div 



V/00V 

2 



(58) 
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As 1//00 satisfies a logarithmic Sobolev inequality with the constant p > 0, we have 
that, for any time t > 0, 

H(y(t,-)\y„) < (2p)- 1 /( V /(f,-)|^)- (59) 
We infer from (|58]l and (59]) that 

Using the Gronwall lemma, we obtain the claimed result. 
Conversely, if 

Vr>0, £{t) < ^(0)exp(-2pj3" 1 f), 



we also have 



and by letting t go to 0, using d58l ), one obtains the logarithmic Sobolev inequality 
I(Y(0,-)\Y«,)>2 P H(y(0,-)\\ir„). □ 



3.2 Effective dynamics 

Consider Q, that solves (O. By a simple Ito computation, we have 

d% (ft ) = ( - W • V| + J3 - 1 A ^ ) (ft ) A + 1 I (a ) (60) 

where is the one-dimensional Brownian motion 

Of course, equation ( f60b is not closed. Following Gyongy fl5l . a simple closing 
procedure is to consider rj, solution to 

dr\ t =b{t,r),)dt + ^j2p- 1 a(t,rj t )dB t , (61) 

where 

b(t,z) =E[(-Vy.V^+j3- 1 4^)(a) |§(ft)=z], (62) 

6r 2 (r,z) =E[|V^| 2 (ft) |§(ft)=z]. (63) 

Note that b and <r depend on f , since these are expected values conditioned on the 
fact that 4 (Q,) = z, and the probability distribution function of ft of course depends 
on t. 



Some remarks on free energy and coarse-graining 



31 



As shown in (15l , this procedure is exact from the point of view of time 
marginals: at any time f, the random variables r/ r and 4(Gr) have the same law. 
This is stated in the following lemma. 

Lemma 8 ([18], Lemma 2.3). The probability distribution function yn of £>{Q t ), 
where Q, satisfies (fJJ, satisfies the Fokker-Planck equation associated to ( 167 1 : 

by/ 6 = d z (-b y/^ +jr 1 a(aV)) • 

The problem with equation d6Tb is that the functions b and a are very compli- 
cated to compute, since they involve the full knowledge of yr. Therefore, one cannot 
consider ( f6Tb as a reasonable closure. A natural simplification is to consider a time- 
independent approximation of the functions b and a. Considering (f62b and (l63l l. we 
introduce (En denoting a mean with respect to the measure jx) 

b(z) = E M [(-W ■ V£ +J3- 1 4^) (Q) | 4(6) = z] 

= J (-W-V^+ZT 1 ^)^, (64) 

and 

a 2 (z) =E P (|V4| 2 (G) \UQ)=z)=J e |V^| 2 rf^, 

where jU^. is defined by (|8}. This simplification especially makes sense if 4 (6«) is a 
slow variable, that is if the characteristic evolution time of 4 (Q t ) is much larger than 
the characteristic time needed by Q, to sample the manifold E-. This is quantified in 
the sequel. 

In the spirit of (loTT i. we next introduce the coarse-grained dynamics 

dlj t =b(l] t )dt+^2^ o(l] t )dB tl lj t=0 = t;(Q ). (65) 

We have proved in [18| that the effective dynamics d65l l is ergodic for the equi- 
librium measure 4 *M> that is exp(— f}A{z))dz. In addition, this measure satisfies 
a detailed balance condition. We have also proved the following error bound, that 
quantifies the "distance" between the probability distribution function of 4(2?) ( at 
any given time t) and that of r\ t . 

Proposition 1 ([181, Proposition 3.1). Assume that 4 is a smooth scalar function 
such that 

forallqeW, < m < |V§(g)| < M < °°, (66) 

and that the conditioned probability measures defined by (0, satisfy a logarith- 
mic Sobolev inequality with a constant p uniform in z: for any probability measure 
V on E z which is absolutely continuous with respect to the measure pL%_, we have 



(67) 
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Let us also assume that the coupling is bounded in the following sense: 

k=\\?z 1 F\\l~«*>, (68) 

where F is the local mean force defined by ( liil) . 

Finally, let us assume that | V£ | is close to a constant on the manifold E z in the 
following sense: 

\Vt;\ 2 -o 2 ot; 



X = 



< oo. (69) 

L°° 



a 2 ot, 

Assume that, at time t = 0, the distribution of the initial conditions of §3§ and ( 1651 ) 
are consistent one with each other: \\f^> (t = 0, •) = (j)(t = 0,-). Then we have the 
following estimate: for any time t > 0, 

M 2 / 2 m 2 B 2 K 2 

- — , A Z H - 

4m z \ p 



E(t) < — 2 A 2 + ) (H( V (0, -)\n) - H( W (t, -)|M)) , (70) 



where E(t) is the relative entropy of the probability distribution function of 
<?{Qt)> where Q t follows 0, with respect to the probability distribution function 
<j) of the solution Tj, to | 



The above proposition thus yields a uniform-in-time bound on the relative en- 
tropy between and <j>. In addition, we also know that the effective dynamics 
is ergodic for exp(— j5A(z))dz, which is the equilibrium measure of £,(Qt), in the 
long-time limit. We thus expect the two probability densities to converge one to 
each other, in the long-time limit. This is indeed the case, as it is shown in lfl8l 
Corollary 3.1]: under some mild assumptions, the L l distance between yn(t,-) and 
(f , ■) vanishes at an exponential rate in the long-time limit. 



3.3 The proof in a simple two-dimensional case 

For the purpose of illustration, we consider in this section an extremely simple 
case: starting from the overdamped dynamics ([3]) in two dimensions (we write 
q = (x,y) € M 2 ), we want to derive an effective dynamics for the coarse-grained 
variable (x,y) = x. Although this case is over-simplified, it turns out that 

the main arguments of our derivation, as well as the proof arguments, can be well 
understood here. 

In that context, the complete dynamics (O reads 
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f dX t = -d,V(X,,Y t )dt+^2j5- l dW t x , 

(71) 



{ dY, = -dyV(X t ,Y t )dt + V2J3 -1 dW?, 
with the initial condition Qq = (Xq,Yq). The manifold E z defined by (0 is 

^ = {(z,y); y^R} 

and the probability measure djj,z, defined by ([H) reads 

d = exp(-pV(z,y))dy = W«(z,y)dy 



e*p(-PV(z,y))dy / V<»{z,y)dy 
Jr 



We focus on the dynamics of £,(X t ,Y t ) ~X t . In that case, the equation ( |6(H > is 
just the first line of ( |7T1 i. which is obviously not closed in X t , since Y t appears. At 
time t, Q t is distributed according to the measure \j/(t,q). Hence, the probability 
distribution function of Y u conditioned to the fact that % (Qt) =X t = x, is given by 

(tv)- V{hXj) 

V / condl'>)7 ~ 



y{t,x,y)dy 



Following Gyongy 1151 . we introduce the function b(t,x) defined by (l62l i. which 
reads in the present context as 



, , d x V(x,y)\j/(t,x,y)dy 
b(t,x) = / [-d x V(x,y)} Wcon d (t,y)dy = -— — j (73) 

/ lift t v v\ 



yf(t,x,y)dy 



and the resulting dynamics (|6TT > reads 



dX, = b(t,X t )dt + \Jl$- x dW t x . (74) 

We now prove Lemma|8]in that specific context and show that, at any time t , the 
probability distribution function of X t is equal to that of ^ (Q t ) = X t . 

Proof (Lemma\8\ case £,(x,y) = x). The probability density function \ir(t,x,y) of 
Qt = (X t ,Y t ) satisfies the Fokker-Planck equation i57Y . 

d,\j/ = div (v/Vy) + j3 _1 4v / 

= d x (\j/d x V) + dy (xi/d,V)+p- l d xx xi/ + l3- l d, y xi/. (75) 

The probability distribution function of £, (Q,) = X t is 



<P(M) = / \jf(t,x,y)dy. 
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Integrating ( 1751 ) with respect to y, we obtain 

d t yfi = djj yr&Vdyj +P~ 1 d xx \^ 



= -d x (ySb)+p- l d xx yS, (76) 



where b(t,x) is given by (1731 . We recognize the Fokker-Planck equation associated 
to the equation d74| i. □ 

As pointed out above, ( l6TT l (i.e. d74l > here) cannot be considered as a reasonable 
closure, since it involves the function b, which is defined using y/(t,x,y) (see (|73ll). 
which in practice is hardly computable. We thus approximate b by the function b 
defined by d64b . which amounts to replacing \j/(t,x,y) in ( |73l by the equilibrium 
measure y/oo(x,y): 

b(x) = -^- . 

/ V»{x,y)dy 
Jr 

In the spirit of d74l i, we thus introduce the effective dynamics 



dX, =b{X t )dt + y/2fi- 1 dW t \ (77) 

We now prove Proposition Q] (error estimator on the effective dynamics), in the 
specific case at hand here. The assumption ( |67] i means that the measure d72t satis- 
fies, for any z, a logarithmic Sobolev inequality with a constant p independent of z. 
The assumption d68l > reads K = \\d xy V\\L a ° < °°, and the assumption (f69t is satisfied 
with A = since V<^ = (l,0) r is a constant vector. 

Proof (Proposition^ case %(x,y) = x). By definition (see the free energy A 
associated to the reaction coordinate % satisfies 

exp(-j3AW)= / Y~(x,y)dy = Z- 1 f exp(-pv(x,y))dy 

hence 

/ d x V(x,y)y«,(x,y)dy 

A'(x) = ± - r =-b(x). (78) 

y/oo(x, y)dy 



The effective dynamics d77l ) thus reads 

d% = -A'(X t )dt + VW<. 

Note that, in this specific context, the effective dynamics is of the form (O (see (TSJ 
Section 2.3] for a comprehensive discussion of the relation between the effective 
dynamics and ©). The probability distribution <j>(t,x) of X t satisfies the Fokker- 
Planck equation associated to the above stochastic differential equation, that reads 
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(t,x)dx. 



We compute, using d79l and d76l i, that 
dE 



dt 



= In 



Jr 



in ^ 



-r 1 / d x 

Jr 

Jr 
Jr 



In 







In 



In 



4> 



yrd x <j> 



b+A 



<p d x 



d T f ln-^J (b+A' 



-p- 1 I(^\(j>) + d x (\n^-\ (b+A 



Using a Young inequality with a parameter a > to be fixed later, we obtain 

. \ 2 







-L -^-1)7(^10) + ?/ i^fA'+fo 



2a 



(80) 



We now observe that, in view of (l73i l and (l78l l. A' and — & are averages of the same 
quantity with respect to different probability measures: 

-b(t,x)= [ d x V{x,y)v^{y)dy and A(x) = [ d x V(x,y) v x (y)dy 
Jr Jr 



with 



V{t,x,y) 
jRV(t,x,y)dy 



and v x (y) 



w°°(x,y) 

J R \j/oo(x,y)dy' 



(81) 
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A'(x)+b(t,x) = / d x V(x,y)v x 2 (y)dy- / d x V(x,y)v[' x (y)dy 
Jr Jr 

= [ (d x V(x, yi )- d x V{x,y z ))kf*(yi,y z ) dyi dy 2 
Jr 2 

for any probability measure k Ux such that 

I k t ' x (y 1 ,y 2 )dy 2 = vi(y 1 ) and f k tx (y u y 2 ) d yi = v' 1 x (y 2 )- 
Jr Jr 

Hence, 

A'(x) + b(t,x) < \\d»V\\ L - [ \yi - y 2 \ ^ (y h y 2 ) d yi dy 2 

Jr 1 

1/2 



< H^VIU- U \yi -nfV'iyun) d yi d y: 



We now optimize on k r ' x . Introducing the Wasserstein distance W(v[ A , V?) between 
v[' x and Vj (see (|55V). we obtain 



A'(x)+b(t,x) ^H^VHi-WCv^vf) 



As recalled above, assumption ( |67l ) means that v 2 satisfies a Logarithmic Sobolev 
inequality. Thus, it also satisfies a Talagrand inequality (see Lemma|7]i, hence 



w(vr,v£) < w^(v^|v|) < -vwiv 



As a consequence, 



A'(x)+A(»,jc) 



< 



I4.V 



Using dSTJ, we obtain 

/ (a' 



< 



< 



\dxyV\\ 



^(t,x)I(v' { x \v x )dx 



d y In 



Y(t,x,y) 



dy 



dx 



Returning to (f80b . and using (l58l l. we thus deduce that 
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We take 2a = j3, so that the first term vanishes, and we are left with 

Integrating this inequality between the times and t, and using that E(0) = 0, we 
obtain 

E(t)< (g(y(f = 0)ly,)-g(y(v)ly-))- 

As recalled above, assumption d68l > reads K" = ||<9 at V||l" < 00 ■ The above bound is 
thus exactly the bound ( |70l i in the present context. □ 



3.4 Numerical results 

In this section, we check the accuracy of the effective dynamics ( |65| l in terms of 
residence times, and also compare this effective dynamics with the coarse-grained 
dynamics © based on the free energy. We perform such comparison on two test- 
cases, and evaluate the influence of the temperature on the results. We also provide 
some analytical explanations for the observed numerical results. 

In the following numerical tests, we focus on residence times. We have indeed 
already underlined that the characteristic behaviour of the dynamics (O is to sam- 
ple a given well of the potential energy, then suddenly hope to another basin, and 
start over. Consequently, an important quantity is the residence time that the system 
spends in the well, before going to another one. 

For all the numerical tests reported in this section, the complete dynamics ([3]) has 
been integrated with the Euler-Maruyama scheme 

X J+1 = Xj - AtVV{Xj) + ^Atp- 1 Gj, 

where, for any j, Gj is a n-dimensional vector, whose coordinates are independent 
and identically distributed (i.i.d.) random variables, distributed according to a nor- 
mal Gaussian law. 

For the simulation of the dynamics (l65l l and ©, we need to have an expres- 
sion for the free energy derivative A' and the functions b and a. These have been 
computed using the algorithm proposed in 171, on a regular grid of some bounded 
interval. Values of the functions for points that do not belong to that grid were ob- 
tained by linear interpolation. We have again used the Euler-Maruyama scheme to 
numerically integrate the dynamics d65b and (O. 
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To compute residence times in a well, we have proceeded as follows (for the sake 
of clarity, we assume in the following that there are only two wells in the test case at 
hand). First, the left and the right wells are defined as the sets {q e R"; | (q) < <^ e h ft } 

and |g € W; %(q) > ^ht} respectively, with ^ r * ht > ^ fv Next, we perform the 
following computations: 

1. we first generated a large number jV of configurations {g,- € R"}i<;<^, dis- 
tributed according to the measure fx restricted to the right well: as a consequence, 

^•)>e g hr 

2. we next ran the dynamics (01 from the initial condition qi, and monitor the 
first time T, at which the system reaches a point q(Ti) in the left well: T, ~ 
inf 

3. from these (fi)i<i</, we computed an average residence time and a confidence 
interval. These figures are the reference figures. 

4. we next consider the initial conditions {^(qi) €E R} 1<;< ,y for the effective dy- 
namics. By construction, these configurations are distributed according to the 
equilibrium measure t, */i (that is exp(— j3A(z)) dz) restricted to the right well. 

5. from these initial conditions, we run the dynamics ((65) or <j6j until the left well 
is reached, and compute, as for the complete description, a residence time and its 
confidence interval. 



3.4.1 A three atom molecule 

Our aim in this section is to show that different reaction coordinates, although sim- 
ilar at first sight, can lead to very different results. As explained in 11181 . the error 
estimate ( TTOb can then help discriminating between these reaction coordinates. 

We consider here a molecule made of three two-dimensional particles, whose 
positions are qA, qB and qc- The potential energy of the system is 

V(q) = ( rAB -l eq f + L (r BC -e eq ) 2 + W 3 (9 ABC ), (82) 

where tab = WlA — Qb\\ is the distance between atoms A and B, £ e q is an equilibrium 
distance, Oabc is the angle formed by the three atoms, and Wi(9) is a three-body 
potential, that we choose here to be a double-well potential: 

W»(0) = hc B ((9- saddle ) 2 - 80 2 ) 2 . 

Wells of W-} are located at 9 — sa ddie ±50. The potential ( l82l l represents stiff bonds 
between particles A and B on the one hand, and B and C on the other hand, with a 
softer term depending on the angle 9abc- To remove rigid body motion invariance, 
we set qB =0 and q& ■ e y = 0. In the following, we work with the parameters e = 
1CT 3 , k e = 208, 4q = 1, ©saddle = ^/2 and 89 — Saddle — 1-187. All dynamics are 
integrated with the time step At = 10~ 3 . 
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We consider two reaction coordinates, that both indicate in which well the system 

is: 

• the angle formed by the three atoms: 

%\ = &ABC- 

In that case, wells are defined by {q g W; % x (q) < and \q g W; (q) > 
with <^ c h ft = 6> saddle - 0.15 and ^ = 0saddle + .15. 

• the square of the distance between A and C: 

& = \\qA-qc\\ 2 - 

In that case, wells are defined by {^elR"; < 4* ft } and {«? g E"; § 2 (?) > ^ 

with = 1.6^ and §* rt = 2.4^. 

Note that there is a region of state space that does not belong to any well. This choice 
allows to circumvent the so-called recrossing problem. 

Remark 8. Note that d82i i reads 

V(q) = (U AB (q) 2 + U BC (q) 2 ) +W 3 (9 ABC ) 

with Uab{q) = r AB — ieq and C/bc(<?) = r BC — 4q- The two first terms in V are much 
stiffer than the last one. We observe that V9abc ' ^Uab = V6Ubc ' ^^bc = 0. Hence, 
the reaction coordinate %\ is orthogonal to the stiff terms of the potential energy, 
in contrast to %%. In view of lfT8l Section 3.2], we hence expect to obtain accurate 
results with §i, in contrast to This is indeed the case, as shown in the sequel of 
this section, o 

We compute the residence time in a given well following the complete descrip- 
tion, and compare it with the result given by a reduced description, based either 
on d65l ) or ||6). Results are gathered in Table Q] for the temperatures J3 _1 = 1 and 
/3 _1 = 0.2. We observe that working with <jji (and either d65l ) or ©) leads to very 
accurate results, independently of the temperature. On the other hand, when the 
reaction coordinate is not orthogonal to the stiff terms of the potential, both coarse- 
grained dynamics turn out to be not accurate. 

Remark 9. In the case at hand here, j|V^i|| 2 = ||V6UbcI| 2 = r Bc- This quantity is 
almost a constant, since the bond length potential is stiff and the temperature is 
small. Hence, along the trajectory, we have that jjV^jl 2 w = 1. This explains 
why, when choosing the reaction coordinate %\, dynamics d65l l and (O give similar 
results, o 

We now study how results depend on temperature. Let us first consider the re- 
action coordinate E,\ = Oabc- Results are shown on Fig. [5] Both coarse-grained dy- 
namics provide extremely accurate results, independently of the temperature. We 
also observe that we can fit the residence time T res according to the relation 
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Temperature 


Reaction 
coordinate 


Reference 
residence time 


Residence time 
using \65\ 


Residence time 
using ((6) 


r' = i 


^1 = & ABC 


0.700 ± 0.011 


0.704 ±0.011 


0.710 ± 0.011 


/H = i 


b = 4c 


0.709 ± 0.015 


0.219 ± 0.004 


2.744 ± 0.056 


J3 ' =0.2 


E,\ = 6 ABC 


5784 ± 101 


5836 ± 100 


5752 ± 101 


J3 ' =0.2 




5833 ± 88 


1373 ± 20 


2135 ± 319 



Table 1 Three-atom molecule: residence times obtained from the complete description (third col- 
umn) and from the reduced descriptions (two last columns), for both reaction coordinates (confi- 
dence intervals have been computed on the basis of jY = 15000 realizations). 



T res «T r ° es exp(sj3) (83) 
with t° = 0.07521 and s = 2.25031. 
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Fig. 5 log jq (residence time) as a function of j8, for the reaction coordinate £i = d^sc- 



By analytical considerations, we now explain why the residence times computed 
from both coarse-grained dynamics © and (l65T l satisfy the relation ([83), with the 
numerical values of s and T r ^. s reported above. 

We first consider the coarse-grained dynamics © driven by the free energy. In 
the case at hand here, it is possible to compute analytically the free energy. Using 
the internal coordinates tab, rgc and 0abc> we indeed infer from (f2]i that the free 
energy A i does not depend on the temperature and satisfies 

A l (e ABC ) = w 3 (e ABC ). 

Thus A i has two global minimizers, separated by a barrier 



AAi = -k e (5e) 4 tt 2.25648. 
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The large deviation theory can be used to understand the behaviour of the dynam- 
ics ([6}, in the low temperature regime. It yields the fact that, when /3 3> 1, residence 
times are given by 



4.D , 



T^ s LU exp(j34Ai) with Tj 



O.LD 



2% 



ft>sp<»w 



(84) 



where ©sp = \J — A'((^sp) is the pulsation at the saddle-point ^sp = ^saddle* and 
ft>W = \M'/(4w) lS me pulsation at the local minimizer £\y = Saddle ±50 (see 
also lfl6l Eqs. (7.9) and (7.10)]). In the present case, we compute that Osp ~ 7.828 



and 0)w ~ 1 1 -07, thus T r 



O.LD 



0.0725, and we find that 



s w A A i and 



^ O.LD 
*res ~ *res 



We thus obtain a good agreement between d83l l and d84l ). as observed on Fig. [5] Note 
that this agreement holds even up to temperature j3 ~ 1 = 1 . 

We now turn to the dynamics d65l ). We pointed out in Remark [9] that dynam- 
ics d65b and (O are identical in the limit of low temperature. The functions b and 
a are plotted for the temperature = 1 on Fig. [6] We observe that, even though 
the temperature is not very small, we already have b w — W 3 ' = —A\ and fj ~ 1 . The 
agreement is even better when the temperature is smaller. This thus explains why 
results given by both coarse-grained dynamics d65l ) and (O can be fitted by the same 
relation d83l ). on the whole range of temperature. 




1.0055 



1.005 - 



1.003 



t — i — i — i — i — i — r 




J I I I I I L 



1 1.2 1.4 1.6 1.8 2 2.2 
Fig. 6 Plot of the functions b and fj, for the reaction coordinate §i = OabC' at the temperature 



We now consider the reaction coordinate = r\ c . Residence times as a func- 
tion of the inverse temperature j3 are shown on Fig. [7] We observe that neither the 
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dynamics © nor the dynamics ((65) provide accurate results. More precisely, the 
reference results, the results given by d65l ) and the results given by (O can be fitted 
by 

^s f ^°ef f exp(,/3), 

T*«r£f exp(*j3), (85) 
^r«^es ree exp(^) (86) 

respectively, with the same parameter s = 2.21 ±0.03 and 

T r °f f « 0.0768, T°' eff « 0.0241 , T Wree « 0.293. 

The dependency with respect to the temperature is thus accurately reproduced by 
both coarse-grained dynamics. The inaccuracy comes from the fact that the prefactor 
T 0,ref j s iii- a pp r0 ximated. 




Again, these numerical observations are in agreement with analytical computa- 
tions based on the large deviation theory. More precisely, we explain in the sequel 
why the residence times computed from both coarse-grained dynamics d65l l and (O 
satisfy ( f83T > and ( |86l l, with the same s, and for the numerical values of s, T°> e and 
T 0,free reported above. 

The functions Ai, b and a are plotted for two different temperatures on Fig. [8] 
Although A2 a priori depends on j3 (as expected), it turns out this dependency be- 
comes quite weak when j3 > 1 . It turns out that we can fit A' 2 by 

A' 2 (§ ) « c 5 {x - 2) 5 + c A {x - 2) 4 + c 3 (x - 2) 3 + c 2 (x - if + Cl (x - 2), 

with a = -16.4433, c 2 = 3.87398, c 3 = 34.2171, c 4 = -6.36938 and c 5 = 
—7.89431. The free energy has thus two local minimizers, ^w,r ~ 2.73 and <^w,i ~ 
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1 .25 and a saddle point, ^sp ~ 2, with 

A 2 (£ SP )«0, A 2 (£ w , r ) « -2.1, A 2 (§ w ,,) « -2.37. 




We introduce the barriers to go from the right well to the left well (r — > 1) and 

vice-versa: 

AAr 1 =A 2 (&p) -A 2 (§ w ,r) and AA^' = A 2 (§ SP ) -A 2 (§ w ,i)- 

In the case of the dynamics © driven by the free energy, and under the assumption 
that the temperature is low enough so that A 2 becomes independent of j3, the large 
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deviation theory can again be used, and yields the fact that residence times are given 
by 

„LD,r-!-l 2ft , M K „LD,l->r 2ft ,„ ■ ,u r , 

T res,free 77 77 ex P(fiAA 2 ), t£ ^ ~ exp(/3AA 2 ), 

where 0)sp, COw.i and «w,r are the pulsations at the saddle-point, the left well and the 
right well, respectively. In the present case, we compute that Osp ~ \f—c\ ~ 4.055, 
c«w,i ~ 5.809 and a)w,r « 4.774. 

The left well is deeper than the right well. Hence, in the low temperature limit, 
the residence time in the left well is much larger than the residence time in the right 
well, and the probability to be in the left well is higher than the probability to be in 
the right well. Hence, 

LD LD.l^r O.LD.l^r /Q , .hri tU O.LD.l^r 2ft 

Cfree « C.ftee ~ Wee exp(/34A 2 ) With T„ s , free = ■ (87) 

®SP f^w.i 

With the parameters that we used, we compute T ~ 0.267, hence 

t ~ AA 1 ^ 1 and T°' free ~ T °' LD1 ^ r 
s ~ ziA 2 ana T res ~ \ es ,fr e e > 

and we obtain a good agreement between (l86*l l and (l87b . 

We now turn to the dynamics d65l l. The functions and <7 plotted on Fig. [8] 
seem to be almost independent of the temperature when j3 > 1. Following lfl8l 
Section 2.3] and ifTTl Sec. 10 and Eq. (89)], we introduce the one-to-one function 

h(t)= a 1 (y) dy and the coordinate £ = h(tz)- We next change of variable in 
Jo 

the effective dynamics (l65T l on the reaction coordinate B, and recast it as 



dC t = -A'(Q)dt + ^2p^dB t , 

where A turns out to be the free energy associated to the reaction coordinate £(q) = 
Mfe (<?))■ The residence time to exit the left well is hence given by 

_LD,l->r 2ft /o i71-»n 

W ~ r r, » exp(/3zU ). 

®SP C0w,l 

In the regime of low temperature, the second term of (fTTb is negligible, and we 
deduce from dTOb that A (/;(<!;)) = A(£,). As a consequence, 

ziA 1 ^ 1 ' = AA 1 ^ 1 ', 5 SP = co SP o-(^sp), Sw.i = fi>w,i cr(| w ,i). 

Hence 
with 
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_0,LD,l-« 27T 



^res,eff 



osp fflwj ct(^sp) c(4w,i) 



We thus recover that the dependency of the residence times with temperature is 
identical between the residence times predicted by the effective dynamics d65l l and 
the residence times predicted by (|6): this dependency is exponential, with the same 
prefactor ZiA 1 ^ 1 '. 

We also compute cr(£ S p) « 3.465 and ct(| w ,i) ~ 2.563, so T^ I->r 0.03. Thus 
the values 1%£ S and qualitatively agree, and we obtain a good agreement 

between (|85) and (|88). 



3.4.2 The butane molecule case 

We now consider a system in higher dimension, namely a butane molecule, in the 
united atom model [24 , 21]. We hence only simulate four particles, whose positions 
are q' for 1 < i < 4. The potential energy reads 

3 

) + Vt>ond-angle(01 ) + ^bond-angle (Ql) + Version (0)> 

i=l 

where 0\ is the angle formed by the three first particles, &> is the angle formed by 
the three last particles, and <p is the dihedral angle, namely the angle between the 
plane on which the three first particles lay and the plane on which the three last 
particles lay, with the convention E {—H, 7r )- We work with 

VbondW = ^{i-ieqf, Vbond-angle(0) = y (9- Q eq f 

and 

Vtorsion(0) = c l(l -COS0) + 2c2(l — COS 2 0)+C3(1 + 3 COS tj) -4cOS 3 0). 

Rigid body motion invariance is removed by setting q 2 = 0, q 1 ■ e z = and q 3 ■ e x = 
q 3 ■ e z = 0. 

In the system of units where the length unit is £q = 1 .53 ■ 10~ 10 m and the energy 
unit is such that k^T = 1 at T = 300 K, the time unit is t = 364 fs, and the numerical 
values of the parameters are l eq = 1, £3 = 208, eq = 1.187, c\ — 1.18, C2 — —0.23, 
and C3 = 2.64. We will work in the sequel with ko = 1000. We set the unit of mass 
such that the mass of each particle is equal to 1 . 

For these values of the parameters c,, the function Vtorsion nas a unique global 
minimum (at 0=0) and two local non-global minima (see Fig. |9). It is hence a 
metastable potential. We choose to work with the dihedral angle as reaction coordi- 
nate: 

§(9) = *. 
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We are interested in the residence time in the main well (around the global mini- 
mizer <j>o = 0) before hoping to any of the two wells around the local minimizers 
(j>±i = ±2n/3. For each minimizer 0o, 0i and the associated well is defined by 
{q; \%(q) - fr\ < £ th }, / = -1,0,1, with ^ = 0.5. 

Remark 10. We observe that 

W stiff -V§=0, 

Where V sm (q) = £- =1 Vbond (\W +l -q'\\) + Vbond-angle(0l)+Hond-angle(02)- In view 

of ifTHl Section 3.2], we hence expect to obtain accurate results with this choice of 
reaction coordinate, as it is indeed the case, o 




Fig. 9 Torsion angle potential Version (0)- 



As in the previous section, we compute reference residence times by integrating 
the complete dynamics, and we then consider both coarse-grained dynamics d65l l 
and ([5). All computations have been done with the time step At = 10~ 3 . Results 
are reported in Table |2] We observe that the effective dynamics (l65l l again yields 
extremely accurate results. The results obtained by the dynamics ©, although qual- 
itatively correct, are less accurate. This conclusion holds for all the temperatures we 
considered. 



Temperature 


Reference 
residence time 


Residence time 
using )65t 


Residence time 
using {6j» 


r' = i 


31.9 ±0.56 


32.0 ± 0.56 


37.1 ± 0.64 


p- 1 =0.67 


493 ± 8 


490 ± 8 


581 ± 9 


=0.5 


7624 ± 113 


7794 ± 115 


9046 ± 133 



Table 2 Butane molecule: residence times obtained from the complete description (second col- 
umn) and from the reduced descriptions (two last columns), at different temperatures (confidence 
intervals have been computed on the basis of ,jY = 13000 realizations). 



As in the previous section, residence times depend on the temperature following 

l^s « T° s exp(>0). 
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For both coarse-grained dynamics, the values found for s and T? es agree with pre- 
dictions based on the large deviation theory. In the case at hand here, it turns 
out that the free energy associated to the reaction coordinate £,(q) = <p is simply 
A (4) = VtorsionOs)' ^ n Fig. [TO] we plot the functions b and a. We observe that they 
are almost independent of the temperature (as soon as j3 > 1), and that a is almost a 
constant. Hence, up to the time rescaling f reS caie = Ot, the effective dynamics reads 
as the dynamics §6^ governed by the free energy. As a = 1 .086 « 1 (see Fig. [Toll, 
the dynamics © yields qualitatively correct results. 




-4 -3 -2 -1 1 2 3 4 -4 -3 -2 -1 1 2 3 4 

Fig. 10 Plot of the functions b and a, for the reaction coordinate ^ = <j>, at different temperatures. 
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